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SEISMIC DATA INTERPRET ATION METHOD 

RACKGRQUNP OF THE INVENTION 

This invention relates generally to seismic data interpretation methods 
and more particularly to a method of identifying subsurface geologic features 
using seismic data. 

The oil industry today recovers only about 35% of the oil it finds. To 
make the next quantum leap in productivity, where ultimate recovery rates 
routinely exceed 50%, the process by which oil and gas reservoirs are 
managed must change. A key tool for these types of improved reservoir 
management activities is the shared earth model, where seismic data is 
integrated with all other available data, such as drilling data, well log data, 
well test data, and production data, to describe the static and dynamic 
properties of the reservoir. 

Most reservoir data is obtained in the immediate vicinity of the well 
bore and measure conditions that are directly associated with only a small 
fraction of the total reservoir volume. Seismic imaging is .typically the only 
method available for identifying the boundary surfaces of hydrocarbon 
reservoirs, including seismic reflector surfaces and faults. This reservoir 
structure information is often crucial for effective field development and 
production planning activities. 

Hydrocarbons typically migrate upward from their source through 
porous subterranean strata until the route is blocked by a layer of 
impermeable rock, and they accumulate beneath this sealing structure. 
Geologists divide traps into two types, structural and stratigraphic. Structural 
traps are formed by tectonic forces after sedimentary rocks have been 
deposited and generally fall into two categories, anticlines and faults. 
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Stratigraphic traps are most often formed at the time the sediments are 
deposited. They fall into three categories: pinchouts, most common in stream 
environments where a channel through a flood plain has been filled with 
permeable sand that was then surrounded by less permeable clays or silts 
when the channel moved; unconformities, where a permeable reservoir rock 
has been truncated and covered by an impermeable layer following a 
depositional period or a time of erosion; and carbonate reefs, fossilized coral 
structures and associated deposits that arose from ancient ocean shelves 
and were overlain by layers of both permeable and impermeable rock. 

In both structural and stratigraphic traps, the containing surfaces of 
hydrocarbon reservoirs invariably consist of interfaces between relatively 
permiable materials and relatively impermiable materials. The bottom 
surfaces of hydrocarbon reservoirs typically comprise hydrocarbon/water 
interfaces located within relatively permiable materials. In many cases, the 
acoustic impedances of materials on opposing sides of these surfaces will be 
significantly different and these boundary surfaces will act as effective 
reflectors of seismic energy. These differences in acoustic impedences 
provide those involved in hydrocarbon exploration and production activities 
with the opportunity to remotely identify hydrocarbon reservoir boundary 
surfaces using reflection seismology. 

Faults are formed when tectonic forces break sedimentary rock, and 
displacements occur along the breaks. A fault may generate a fluid channel 
through a sealing layer, or seal a permeable layer. 

A conventional approach to identifying both seismic reflectors and fault 
surfaces using seismic data is to view two dimensional cross sections of the 
seismic data and to manually identify either points that appear to lie on a 
common seismic reflector or points at which the primary geologic layers 
appear to be displaced (i.e. a fault). While it may be a relatively 
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straightforward problem to manually interpret seismic reflectors or faults on a 
single two dimensional seismic section, it becomes an extremely tedious and 
difficult task to manually identify these geologic features as three dimensional 
(i.e. planar) surfaces because the number of points that must be identified 
goes up exponentially. 

This is particularly true when identifying seismic reflector surfaces in 
stratigraphic traps. Exploration seismologists find structural traps far easier to 
identify than stratigraphic traps in both 2D and 3D seismic data because 
structural traps are seen as highly dipping reflectors and discontinuities in 
otherwise smooth reflections. For years, seismic data acquisition and 
processing techniques have been tailored to accentuate these features, 
allowing interpreters to concentrate their efforts on faults and anticlines rather 
than their more subtle stratigraphic counterparts. 

As a majority of the major hydrocarbon production areas have already 
been explored for easily identified structural traps, an increasing emphasis is 
being placed on locating and evaluating complicated structural and 
stratigraphic traps. A conventional approach to identifying stratigraphic traps 
is to perform a manual interpretation of all the primary layers, analyze their 
termination points, and then to manually classify the type of trap based on the 
interpreter's training and experience. This is often an extremely tedious and 
difficult task which requires substantial geological experience to allow 
prospective stratigraphic traps to be properly identified and analyzed. 

The complexity of identifying these types of traps using seismic data is 
compounded by changes in the type of seismic data being acquired. The 
increasing use of three-dimensional ("3D") acquisition techniques and high- 
resolution surveys have dramatically increased the quantity of the seismic 
data available to be interpreted. At the same time, stable to falling crude oil 
prices and competitive pressures have required companies to more efficiently 
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utilise their exploration and production expenditures and to reduce the time it 
takes to reach decisions regarding the development of new fields and 
prospects. 

Fortunately, modern computer systems offer particularly powerful tools 
for processing and aiding in the interpretation of seismic data. While the 
processing capabilities of computer systems have been increasing 
exponentially, the cost of these systems have simultaneously fallen 
dramatically. Modern computer hardware and software systems are being 
used successfully to process seismic data and to produce clearer and more 
detailed maps of the subsurface, which often eases the process of manual 
seismic data interpretation. 

The continuing reliance on the manual interpretation of seismic data, 
however, both adds significant costs to and introduces substantial uncertainty 
into the process of seismic data interpretation. Even the most highly 
regarded seismic data interpreter is capable of making errors in judgement or 
suffering from lapses in concentration. 

An even more pervasive problem with manual seismic data 
interpretation is that a human being is efficient at working on only one 
problem at a time and is able to only work for limited periods of time before 
the quality of the work performed begins to rapidly degrade. If an automated 
computer-based seismic data interpretation system is used, however, large 
quantities of seismic data can be subdivided and the task of interpreting each 
portion of the seismic data can be delegated to a separate computer system 
which can be run 24 hours a day, 7 days a week until the process is 
complete. This not only greatly speeds up the seismic data interpretation 
process, but it also eliminates the variation in results that would be inherent if 
more than one person was involved with interpreting the seismic data. 
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Prior art techniques for identifying subsurface geologic features using 
seismic data have typically been limited to identifying features in two- 
dimensional cross-sections. Some of these systems have required digitized 
horizons or operator selected points on the subsurface feature, such as U.S. 
Patent No. 5,229,976, issued to Boyd et al. on July 20, 1993. Other systems 
have used positive/negative amplitude cross-over points and edge following 
techniques to track horizons, such as U.S. Patent No. 5,148,494 issued to 
Keskes on September 15, 1992. 

These types of prior art methods have suffered from three primary 
problems. First, these methods have only identified a linear cross-sectional 
view of the subsurface feature (instead of providing the planar surface 
information desired for shared earth reservoir modelling). Second, these 
identification techniques only utilised two-dimensional input data {limiting the 
accuracy and effectiveness of these prior art methods). Third, the input 
seismic data was required to exhibit strong contrasts or the method would fail 
to properly identify the desired geologic feature. 

It is therefore desirable to implement an improved method of 
interpreting seismic data that overcomes one or more of the problems 
exhibited by these prior art methods. 

An object of the present invention is to provide an improved method of 
interpreting seismic data. 

An advantage of the present invention is that three-dimensional 
geologic features may be identified and the identification process may utilize 
seismic data from three-dimensional neighborhoods of cells to identify the 
geologic feature. 
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Another advantage of the present invention is that it is capable of 
identifying subtle geologic features. 

An additional advantage of the present invention is that the selected 
cells may be clustered to distinguish between different geologic features. 

A further advantage of the present invention is that the identified 
geologic feature may be represented as a plane or other type of polynomial 
surface using a vectorizing procedure. 

Another advantage of the present invention is that the geologic 
features identified may comprise hydrocarbon reservoir boundaries for use in 
three-dimensional shared earth reservoir modeling. 

An additional advantage of the present invention is that the method 
may be substantially automated and may be used to identify subsurface 
geologic features with little or no manual operator input. 



SUMMARY OF THF INVENTION 

The present invention involves a method of identifying a subsurface 
geologic feature using seismic data. The method includes the steps of: 
subdividing a subsurface area into a plurality of discrete cells; associating the 
cells with seismic data that image subsurface regions represented by the 
cells; determining structural characteristics for a plurality of the cells; and 
grouping together nearby cells based on the structural characteristics to 
identify the subsurface geologic feature. The grouped cells may be suitably 
clustered and/or vectorized. The interpretation results may be suitably 
represented and/or displayed. The invention and its benefits will be better 
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understood with reference to the detailed description below and the 
accompanying figures. 



RR1EF DESCRIPTION OF T HE DRAWINGS 

Figure 1 is a process flow chart showing steps associated with the 
inventive method; 

Figure 2 is a schematic plan view of a seismic survey vessel and 
associated equipment acquiring seismic data; 

Figure 3 is a schematic diagram of a seismic data interpretation 
workstation; 

Figure 4 is a cross-sectional view through a seismic data volume with 
highlighted lines showing faults intersecting the section; and 

Figure 5 is a three dimensional perspective view of faults shown in 
Figure 4. 

Figure 6 is a cross-sectional view through a seismic data volume with 
highlighted lines showing seismic reflectors intersecting the section; 

Figure 7 is a three dimensional perspective view of seismic reflectors 
shown in Figure 6; 

Figure 8 is a cross-sectional view through a termination attribute cube 
in accordance with an embodiment of the present invention; and 
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Figure 9 is a cross-sectional view through an attribute cube that 
highlights chaotic reflector responses in accordance with an embodiment of 
the present invention. 



DFTAILED DESCRIPTION OF THE INVENTION 

Figure 1 is a process flow chart showing various steps associated with 
the present method of identifying subsurface geologic features using seismic 
data. In the Acquire Seismic Data step 10, seismic data is obtained from a 
subsurface area. This seismic data is then processed in the Process Seismic 
Data step 12 to improve the signal to noise ratio of the data and to produce a 
more accurate representation of the subsurface geology. The subsurface 
area of interest is subdivided into a plurality of discrete cells in the Subdivide 
Area step 14, with each cell representing a different region within the 
subsurface area. The cells are then associated (i.e. "loaded") with seismic 
data that image the subsurface regions represented by the cells in the Load 
Seismic Data step 16. 

In the Determine Structural Characteristics step 18, structural 
characteristics are determined for a number of the cells. Nearby cells 
identifying the subsurface geologic feature of interest are grouped together in 
the Group Cells step 20. If desired, the grouped cells can then be clustered, 
represented mathematically, such as by a planar surface, and/or displayed in 
the Represent/Display Results step 22. The individual steps associated with 
the present method will now be discussed in detail. 

The activities typically performed in the Acquire Seismic Data step 10 
are discussed with reference to the ship and associated equipment shown in 
Figure 2. In Figure 2, a seismic survey vessel 30 is shown towing three 
seismic sources 32, such as airguns or arrays of airguns. The survey vessel 



WO 99/64896 



PCT/IB99/01040 



9 

30 is also shown towing a pair of seismic streamers 34 that contain large 
numbers of seismic sensors, such as hydrophones. Acoustic pulses 
produced by the seismic sources 32 are transmitted through the water and 
into the geologic subsurface. When the spreading acoustic pulses reach 
subsurface interfaces where the acoustic impedances of materials on 
opposing sides of the interface change, a portion of the acoustic energy is 
reflected and returned toward the surface (i.e. the subsurface interfaces act 
as seismic reflectors). Seismic sensors within the seismic streamers 34 
receive this reflected acoustic energy. The acoustic energy received by the 
seismic sensors is measured (typically at a regular time-based sampling 
interval, such as every 2 milliseconds), digitized and transmitted to the survey 
vessel 30 where the measurement values are recorded. These types of 
sensed, digitized, and recorded acoustic energy-based measurement values 
typically comprise the type of seismic data used in connection with the 
present method. 

Preferably, the seismic data is acquired using a seismic data 
acquisition system that individually digitizes and processes the acoustic 
energy received by each seismic sensor, such as the type of system 
described in PCT International Application No. PCT/GB97/02544, 
International Publication No. WO 98/14800. This type of "single sensor" 
seismic data acquisition equipment is capable of acquiring seismic data that 
provides significantly clearer images of geologic features than conventional 
coarse spatial sampling "hardwired array" types of seismic data acquisition 
systems. 

The seismic data will also preferably be acquired using a high 
resolution, high fold seismic data acquisition geometry. Using such an 
acquisition geometry, the seismic data may be acquired at 12.5m by 12.5m, 
or even 6.25m by 6.25m, bin sizes and the number of independent 
measurements obtained at each bin (the "fold" of the seismic survey) may be 
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48 or higher. Seismic data acquired under these conditions will typically allow 
more consistent and more accurate identifications of geologic features to be 
made using the present method. 

Time lapsed or '^D" seismic data may also be used in connection with 
the inventive method. Acquiring seismic images of hydrocarbon reservoirs at 
different times can provide important information for reservoir characterization 
and monitoring activities. When used with time lapsed seismic data, the 
present method can be used to identify fluid/fluid interface surfaces in the 
hydrocarbon reservoirs (oil/water and gas/oil interface surfaces) and to 
monitor the movements of these surfaces due to hydrocarbon production 
activities. This type of monitoring can provide invaluable information to a 
reservoir engineer contemplating alternative production enhancement and 
remedial well activities. 

Other types of seismic data (in addition to conventional 
pressure/pressure transmission mode seismic data) may also be used with 
the present method to identify subsurface geologic features, such as 
pressure/shear transmission mode seismic data, shear/shear transmission 
mode seismic data and multi-component seismic data. 

The vessel and equipment configuration shown in Figure 2 merely 
illustrates one of a vast number of possible seismic data acquisition systems 
that may be used in connection with the inventive method in marine, transition 
zone, or land environments. The applicability and scope of the present 
method is in no way limited or restricted to the use of one particular type of 
seismic data or seismic data acquired by one particular type of seismic data 
acquisition system. 

After being acquired, the seismic data is processed in the Process 
Seismic Data step 12. This processing typically improves the signal to noise 
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ratio of the data and transforms the seismic data into a more accurate 
representation of the subsurface geology. Typical processing steps may 
include migration (such as pre-stack depth migration), stacking, and/or 
filtering (such as K-F or tau-p domain filtering). The use of decomposed and 
reconstructed seismic data, such as the type of decomposed and 
reconstructed seismic data disclosed in PCT International Application Number 
PCT/IB98/00209, International Publication Number WO 98/37437, 
incorporated herein by reference, may be particularly useful when identifying 
extremely subtle subsurface geologic features. The seismic data may also be 
converted from the time domain (i.e. where the sampling interval of the 
seismic data is a time interval, such as 2 milliseconds) to the depth domain 
(i.e. where the sampling interval of the seismic data is a distance interval, 
such as 1 meter). Interpreting seismic data after it has been converted to the 
depth domain allows the identified features to be more readily integrated into 
a shared earth reservoir model. 

Figure 3 depicts components of a seismic data interpretation 
workstation computer, including one or more portable storage devices 40, a 
hard storage device 42, a CPU 44, one or more operator input devices (such 
as a keyboard 46), and one or more output devices (such as a display 48). 
The workstation may comprise, for instance, a general purpose Silicon 
Graphics lndigo2 workstation computer. The portable storage device 40 may 
comprise magnetic or optical recording media, such as floppy disks, CD- 
ROMs, or magnetic tapes. The portable storage devices 40 may contain 
seismic data or computer software instructions that are copied to the hard 
storage device 42 and that allow the workstation to perform one or more of 
the steps or processes associated with the present method. 

The method steps that follow the Process Seismic Data step 12 are 
preferably performed using the GeoFrame software package available from 
GeoQuest, a division of Schlumberger Technology Corporation, Houston, 
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Texas, USA. This type of software package is particularly suited for 
subdividing the subsurface area of interest into a plurality of discrete cells (the 
Subdivide Area step 14) and associating the cells with portions of the seismic 
data that image subsurface regions represented by the cells (the Load 
Seismic Data step 16). 

In contrast to some prior art methods, the cells in the present method 
preferably contain seismic data associated with a single subsurface point (i.e. 
the seismic data images a subsurface area nearer the centerpoint of this 
particular cell than the centerpoint of any other particular cell), not seismic 
data traces (seismic data associated with a suite or group of spaced apart 
subsurface points). In many cases, the single seismic data value in the cell 
will be a zero-offset pressure/pressure transmission mode seismic amplitude 
value associated with the cell centerpoint. In other cases, more than one 
seismic data value may be associated with a cell, such as a pressure/shear 
transmission mode value in addition to the pressure/pressure transmission 
mode value discussed above, or one or more non-zero offset value 
associated with the imaging point. 

Next, structural characteristics are determined for at least some of the 
cells in the Determine Structural Characteristics step 18. A choice must 
typically be made at this point whether to perform a true three-dimensional 
("3D") identification of the subsurface geologic features or to identify the 
features using more simplistic two-dimensional ("2D M ) methods. Subsurface 
feature identification using 3D methods provides better results, but is more 
computationally intensive. A faster way to accomplish similar results is to 
identify the subsurface features in 2D slices of the data and then to aggregate 
the results. This substantially reduces the computational requirements of the 
method, but the output of the process may not be as accurate or reliable. 
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If three-dimensional calculations are made, the structural 
characteristics determined may comprise or include local seismic reflector dip 
and azimuth values. If two-dimensional calculations are made, the structural 
characteristics may comprise or include seismic reflector local plunge 
estimates (i.e. estimates of the dip angles of the seismic reflectors from the 
particular cross-sectional reference viewpoint of the section). The confidence 
or likelihood that the values of these determined structural characteristics are 
accurate can also be quantified during this process. 



3D Structural Characteristics Determination Methods 



Various three-dimensional methods for determining the local 
dominating dip and azimuth of the seismic reflector will now be described. In 
the following, if a signal, xft,f 2 ,y, has a fixed dip and azimuth, then the loci of 
the signal's iso-value surfaces, i.e. the points where the signal has constant 
value, will be parallel planes. Furthermore, the gradient of the signal, V* , will 
be perpendicular to these parallel planes. 

In a real world dataset, the loci of the iso-value surfaces are not likely 
to be exact planes. Hence, we attempt to estimate the dominating dip and 
azimuth. Given a set of 3-dimensional spatial gradient vectors. 



Vx(*,,/ 2 ,/ 3 ) = 



&('i,' 2 >'3) 



5 *(/ p / 2 ,/ 3 ) 



*3 



the dominating dip and azimuth may be determined by principal component 
analysis. In principal component analysis, the covariance matrix, 
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of the gradient vectors is constructed, and the dominating dip and azimuth 
selected to be given by the principal eigenvector of this matrix. The C»s of the 
covariance matrix are defined by 



(dx ' 



where 



»7* 



Having only one dip and azimuth estimate for an entire signal may be 
of limited practical Utility. To increase the usefulness of this estimate, the 
dominating local dip and azimuth is needed. A localized estimate is obtained 
by replacing the global covariance estimate by a windowed local estimate, 



C»('l>'2>' 3 ) = 



'I .'2 .'3 



ax(/ l 4-r t ,/ 2 +r 2 J 3 +r 3 ) Y dx{t x + r„/ 2 + r 2 ,/ 3 + r 3 ) ' 

5; H 



The window function, w(s) = w(t 1t tM wHl typically be a low-pass filter. 
An implication of the windowing function is a smoothing of the dip and 
azimuth estimate. 



In summary, the local dip and azimuth estimation is done by: 
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1 . Orientation selective filtering, typically by gradient estimation; 

2. Windowed covariance estimation of the gradient vectors; and 

3. Principal component analysis. 

There are several possible implementations of the gradient estimation 
and window function. Below is an exemplification of particularly useful choices 
of these operators. 

Orientation selective filtering: For discrete data, the gradient is a 
discrete estimate. By nature, gradient estimation is rather noise sensitive. In 
order to allow some noise removal, the gradient estimation is done by filtering 
with a derivative of Gaussian ("DoG") filter. If discretization effects are 
ignored, this is equivalent to smoothing the data by a Gaussian low-pass filter 
and then differentiation. The smoothing removes noise, and the size of the 
Gaussian low-pass filter determines the degree of noise removal. 

The unit pulse response of a multi-dimensional DoG filter is separable, 
and has the equation 




in the differentiated dimension and 
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in non-differentiated dimensions, where C K and C k are scalar factors. A 3D 
DoG operation consists of applying h r (k) once and h x {l) twice, to different 
dimensions. For example applying h K (k) vertically and h x (l) to the two 
horizontal dimensions gives the vertical gradient component. The result of 
applying all DoG filters is a gradient vector, Vx(/,,/ 2 ,/ 3 ) . with one partial 
derivative component for each dimension. Other types of orientation selective 
filters may be used to produce orientation estimates that may be used instead 
of or in addition to the gradient estimates described above, including band- 
pass or high-pass filters having orientation selective properties. 

Windowed covariance estimation: At strong edges in the data, the dip 
and azimuth of the gradient, Vjc(/,,/ 2 ,/ 3 ) , will be consistent, but in regions of 
the data with weaker signal, the gradient will have less consistent dip and 
azimuth. This problem is primarily dealt with by the window function in the 
covariance estimate, which essentially smoothes the dip and azimuth 
estimates. 

A practical implementation of this operation is to estimate the volume 



'ac(r,,/ 2 ,r 3 ) Y dx(/,,/ 2 ,f 3 ) 



dt, 



~1j 



and then smooth it with a Gaussian low-pass filter, 



*(/„/,./,) = • 



exp 



( 

1 






2\ 


2 




a 7 <7 3 ; 





The samples of the smoothed volume are then used as the estimates of the 
CfS. This implies the computation of 3x3=9 volumes. Note that the Gaussian 
filter is separable, which speeds up the computation significantly. 
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Principal component analysis: Now we typically have one covariance 
matrix C for each position in the data volume. The dominating local dip and 
azimuth is estimated by calculating the principal eigenvector of C. Detailed 
information regarding suitable principal component analysis implementation 
methods may be found in R.A. Johnson and D.W. Wichern. Applied 
Multivariate Statistical Analysis. 3rd ed. Prentice-Hall International, Inc., 
Englewood Cliffs, NJ, 1992, incorporated herein by reference. 

Estimate reliability measure: There will be three eigenvectors, v,., of the 
matrix C, each of them associated with one eigenvalue, A f . The larger v., the 
better describes the dip and azimuth. The larger the difference between 
the dominating k k and the two other A/s, the more reliable the dip and 
azimuth estimate is. Assuming without loss of generality that A, > X 2 > A 3 , 
one type of calculated value that quantifies the estimates' reliability is: 

Other types of estimate reliability measures may also be utilised in connection 
with the inventive method. 

Parameter tuning: Empirically, it has been determined that useful 
values for a, may be in the range of 0.5 to 3.0 and for cr 2 may be in the 
range of 1 .5 to 6.0. Too large coefficients may blur the dip and azimuth 
estimates too much, while too small values may make the estimates too 
noisy. The parameter values of a, = 0.5 and for a 2 = 3.0 have been 
empirically determined to yield good results in some applications for several 
data sets, including the data set shown in Figure 8. Other applications may 
require different values. 
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2D Structural Characteristics Determination Methods 

The dip and azimuth estimation approach described thus far for three- 
dimensional data may also be specialised for two-dimensional signals. It can 
be shown that the two-dimensional specialisation of the above procedure 
simplifies to the approach presented below. 

The gradient vector is estimated similarly as in the three-dimensional 
case, i.e. by orientation selective filtering, such as by DoG filtering. The angle 
of this vector will have a good correspondence with the angle of plunge in the 
image at strong edges. However, in regions with low energy (like in between 
seismic reflectors), the amplitude of V* will be low and the angle estimate 
less reliable. Hence, some smoothing of the gradient angle, weighted by the 
gradient magnitude, is desired. Note that this smoothing corresponds to the 
windowing in 3D. 

However, adjacent vectors at angles separated by approximately n 
(i.e. pointing in opposite directions) with similar magnitudes will tend to cancel 
each other with such a smoothing. To alleviate this problem, the vectors are 
transformed to a space where vectors separated by an angle of n are 
identical. This may be done by representing V* by complex numbers, with 
e.g. the horizontal component corresponding to the real part and the vertical 
to the imaginary. Then the arguments of the squared complex numbers are 
smoothed, weighted by their amplitudes. To illustrate the effect of this 
process, consider the two complex numbers 
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and 

Squaring yields 

z 2 x =2 2 e J24> ' 50 

and 

z 2 2 = 2 2 e j2{05 * ;r) = 2 2 e JiU2n) =2 2 e J(2050) =z 2 i . 

Hence, z 2 \ and z 2 2 are equal. 

Finally, the estimate of the angle of plunge is obtained as the angle of 
the square root of the elements from the process above. Note that with 
seismic data, it is meaningless to consider angles separated by n as different 
(e.g. a horizontal reflector has no direction left-to-right or right-to-left). 

As in 3D, several possible realizations of the smoothing and derivative 
filters exist. Sets of particularly good candidates are Gaussian and DoG 
filters. As we have seen, these filters are separable, lending themselves to 
particularly efficient implementations. 

With these choices, the 2D angle of plunge estimator may be 
summarized as: 

1 . Orientation selective filtering, typically gradient estimation by 
filtering with a vertical and a horizontal partially differentiated Gaussian filter. 
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2. Combine these two partial derivatives into a vector. 

3. Transform the resulting vector elements by squaring the 
complex representation. 

A. Weighted Gaussian smoothing (accomplished by smoothing the 
complex image). 

5. Square root of the complex result. 

6. The angles of these complex numbers correspond to the local 
angle of plunge. 

This approach corresponds to the three-dimensional approach in that 
steps 1 and 2 correspond to the gradient estimation process and steps 3, 4, 
5, and 6 correspond to the windowed covariance estimation and principal 
component selection processes. 

Other methods for determining similar structural characteristics are 
known and may be used in connection with the present method, including 
those described in T.P.H. Steeghs, Local Power Spectra and Seismic 
Interpretation, Ph.D. dissertation, Delft University of Technology, Delft, the 
Netherlands, Sept. 1997 and M.S. Bahorich and S. Farmer, "3-D seismic for 
faults and stratigraphic features: The coherence cube" The Leading Edge, 
vol. 14, pp. 1053-1058, Oct. 1995, both incorporated herein by reference. 

Fault Structural Characteristics Determination Methods 

The types of structural characteristics described thus far have been 
associated with seismic reflectors, typically the interface surfaces between 
different geologic sedimentary layers. If the subsurface geologic feature to be 
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identified by the present method is a fault, these types of structural 
characteristics can be used to calculate structural characteristics associated 
with subsurface fault features. These fault structural characteristics 
determination methods are based on the dip and azimuth estimation 
technique described above. 

If a fault intersects a strong seismic reflector, the reflector typically 
becomes discontinuous and is characterised by an abrupt change in the 
seismic signal. Derivative filters will enhance abrupt changes, while 
suppressing smooth regions. However, the reflection layers are also 
intrinsically abrupt signal changes, and will also be enhanced. Hence it is 
necessary to perform the differentiation of the seismic signal along the 
reflection layers. In general, it can not be assumed that the layers are 
horizontal. The dip and azimuth of the layers therefore need to be estimated 
first. 

A separable differentiation may be performed, where the components 
of the derivative are combined according to the dip and azimuth estimates. 
For example, if the angle of plunge of the seismic at a particular position in a 
2D section is 6 and the vertical and horizontal partial derivatives are 

dx{n^n 2 ) 

and 



respectively, a structural characteristic that quantifies the likelihood that a 
fault feature passes through a particular cell is 
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dn x dn 2 

This type of structural characteristic is hereinafter referred to as a "fault 
characteristic". 

In 3D, the local dip and an azimuth value represent a plane. 
Consequently, the fault characteristic may be interpreted as the projection of 
the vector with the three partial derivatives, 

dx(n l9 n 29 ni) 
dn } 

Vjc= dx{n l9 n 79 n y ) 
dn 2 
dx(n ]9 n 2 ,n 3 ) 
dn 2 

onto the dip/azimuth plane. In smooth regions, the projected vector will have 
a small magnitude, whereas it will be larger near abrupt signal changes (i.e. 
near a fault feature). 

Fault Cell Grouping/Thinning 

The fault structural characteristics values will have high values when it 
appears that a fault passes through a particular cell and low values when it 
appears that a fault does not pass through a cell. In its most simplistic form, 
fault cell grouping may comprise grouping or associating cells having 
relatively high fault structural characteristic values with nearby cells also 
having relatively high fault structural characteristic values. This step is shown 
generally in Figure 1 as the Group Cells step 20. Three primary problems, 
however, may complicate the fault feature identification process: 
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1. The faults may have discontinuous structural characteristic 
values. That is, the fault structural characteristic values may occur as a series 
of blobs instead of as a continuous region. This appears to be primarily due to 
the lack of appropriate seismic signals in areas between seismic reflectors; 

2. The noise behind the signal may distort the determined seismic 
characteristics; and 

3. The seismic response normal to the fault may be too weak for 
easy fault feature identification. 

An experienced human interpreter may be able to pick the centre-line 
of the determined fault structural characteristic values, but for a computer to 
make a robust pick, a "thinning" or reduction procedure should typically be 
applied. 

Smoothing the fault attribute will reduce problems associated with 
noise and lack of continuity. On the other hand, unless the smoothing filter is 
well aligned with the dip and azimuth of the fault, it would increase the 
problem with blurred responses. However, by smoothing, it can be expected 
that the peak of the smoothed attribute will be close to the center of the fault 
response. 

To utilize this, the neighborhood of the fault is traced perpendicular or 
nearly perpendicular to the fault. Assume that we have dip and azimuth 
estimates in the vicinity of a fault that are perpendicular or nearly 
perpendicular to the fault and "flow lines" are lines following plunge angles or 
local dip and azimuth estimates along its path. 



WO 99/64896 



PCT/IB99/01040 



24 

For 3D attributes, a new dip and azimuth estimate is therefore 
necessary, estimating the orientation perpendicular to the fault. One way to 
make this estimate is to use the dip and azimuth of the gradient vector after 
projection onto the orientation plane. This vector will generally be semi- 
orthogonal to the faults. In order to have reliable dip and azimuth estimates, 
these estimates should also be smoothed. 

Consider, for instance, a trace of fault characteristic values extending n 
samples to each side of the current position. During the thinning process, if 
the value at the current position is the largest value in the group, it is retained 
and otherwise it is discarded. As a result, the fault characteristic response is 
thinned and only cells with the largest amplitudes are retained. This type of 
process may also be used to "thin" or "depopulate" other types of structural 
characteristics discussed below. 

In summary, the process of computing a 3D set of fault structural 
characteristics using this embodiment of the present method may be 
described as follows: 

1 . Read a 3D seismic sub-volume from a cellular seismic 
database. 

2. At each cell, estimate the local dip and azimuth of the seismic 
reflectors, corresponding to a local plane. This plane is represented by its 
normal vector. The dip and azimuth estimate includes the selection of two 
Gaussian filter parameters. 

3. Compute the gradient (more precisely the derivative of 
Gaussian or "DoG") in the seismic volume. The width of the Gaussian is a 
parameter of this step. To have a sharp response, this parameter should be 
quite low. 
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4. Project this gradient vector onto the local plane. This projection 
is now an initial fault structural characteristic, where the projection removes 
the derivative component caused by the reflector. The only components left 
are the components due to faults and other discontinuities in the signal. 

5. Smooth the dip and azimuth of the projected derivative as 
described. It is not very critical to have a dip and azimuth estimate that 
provides a very accurate fault plane orientation. The use of this parameter is 
for tracing a flow line across the fault, and deviations from the "normar or 
"perpendicular" direction can be tolerated. It is, however, important that the 
dip and azimuth is consistently non-parallel to the fault. Consequently, the 
smoothing filter should be quite large. 

6. Compute a refined fault structural characteristic by Gaussian 
smoothing of the magnitude of the result from step 4. The smoothing 
removes noise and "connects" discontinuous fault responses. The size of this 
smoothing filter, cr 5 , is probably the most critical parameter of this approach 
and the parameter can be regarded a scale-space parameter. Consequently, 
a small a 5 corresponds to a fine scale, and subtle fault responses and fine 
detail is enhanced. However, with a small a 5 , the noise suppression is low, 
and a lot of minor "artificial" faults are identified. A large cr 5 , on the other 
hand, suppresses more noise and only enhances the faults with major 
discontinuity responses. 

7. Compute the displacement corresponding to the dip and 
azimuth above. 



WO R9/64896 



PCT/IB99/01040 



26 

8. Trace flow lines beginning at each cell. If the beginning cell of 
the flow line has the maximum structural characteristic value, this value is 
retained. Otherwise, it is marked as blank. 

9. Write the resulting sub-volume to the seismic database. 

Other types of subsurface geologic features may also be identified by 
the inventive method, such as seismic reflectors. 

Common Reflector Cell Selection/Thinning 

In its most simplistic form, cells associated with a common reflector 
may be selected merely by selecting cells having certain seismic reflector 
inclination values or estimate reliability measures. Groups of nearby cells in a 
time to depth converted seismic data set having zero or nearly-zero seismic 
reflector inclinations (i.e. lying on a horizontal or nearly horizontal surface) 
may, for instance, be selected. These types of seismic reflectors may 
comprise hydrocarbon/water interfaces and may be direct predictors of 
hydrocarbon deposits. These types of "flat spots" may be readily extracted 
from seismic data using the inventive method under certain combinations of 
subsurface geology and fluid fill conditions. 

Nearby cells having relatively high estimate reliability measures can 
also be selected. Cells that have relatively high estimate reliability measures 
are typically located on the edges of locally dominant seismic reflectors. 

In some cases, the selected seismic reflector cells will need to be 
thinned. An experienced human interpreter may be able to pick the "center of 
gravity" of a group of selected cells, but for a computer to make a robust pick, 
a "thinning" or reduction procedure may need to be applied. The reflector dip 
and azimuth estimates may be used, for instance, to calculate a vector 
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largest value in the group, it is retained and otherwise it is discarded. 
Alternatively, the center-most cell in the group may be retained. This type of 
process may also be used to "thin" or "depopulate" selected termination cells, 
as discussed below. 

Seismic Reflector Extraction by Flow Surface Generation 

As described above, in some cases seismic reflector positions may be 
extracted from seismic data directly on the basis of particular seismic reflector 
inclination values or estimate reliability measures. In other cases, the Group 
Cells step 20 will involve the identification of flow surfaces. Starting from a 
seed position, a "flow surface" is the surface through that position having the 
angle of the local dip and azimuth estimates (perpendicular to the gradient 
angle) in its path. Flow surfaces may be 2D surfaces in a 3D volume, 1D 
lines in a 3D volume, or 1D lines on a 2D cross section. In the case of lines, 
the flow surface may be denoted a flow line. Flow lines on a cross section 
may be generated by using a 2D angle of plunge estimate, or by projecting 
the dip and azimuth onto the section. A flow line is a straight-forward 
specialization of a flow surface and will not be treated separately. A flow 
surface may be represented as a time (or depth) interpretation, i.e. as an 
interpretation grid, or as a set of surface patches. 

The inventive method is particularly suited for extracting subtle seismic 
reflectors that are difficult or impossible to extract using conventional 
methods. Because the inventive method calculates a local seismic reflector 
inclination for each of the cells examined, a flow surface associated with each 
cell in the database may be generated, regardless of seed cell and the 
relative magnitude of the seismic response from that cell. Many prior art 
methods will fail to identify any surface whatsoever if the magnitude of the 
contrast in the seismic data is small. 
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Seed Cell Selection 

In some cases, seed cells may be manually input because they 
represent known points on a subsurface geologic feature, such as known 
reservoir boundary intersection points. Cells associated with beginning and 
end points of an oil column could be determined from a well log and input in 
the present method as seed cells to produce a detailed image of the reservoir 
cap rock and the oil/water interface, for instance. 

In other cases, it may be desirable to identify a large number of 
seismic reflectors within a data set. Two seed cell layout schemes have 
proven useful in this type of application of the inventive method. 

1 . Seed cells are located on a set of traces with several seeds per 

trace. 

2. Having some scanning order of the samples, additional seed 
cells are located at positions where the distance to the closest previously 
identified geologic feature is above some threshold. 

It should be remembered that the selection of one or more seed cells 
merely allows the flow surface cell selection process a place to start. The 
seed cells can therefore also be thought of as beginning locations or areas in 
the cell grouping process. 
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Momentum 

The concept of "momentum" can also be used during the flow surface 
generation process. When using this momentum component, the local dip 
and azimuth estimates are smoothed over neighboring estimates on the same 
surface. Hence, if the angular estimates would make a flow surface bend 
using the previous approach, now the flow surface is forced to continue 
straighten 

This flow surface generation process may be terminated when the flow 
surface meets a previously generated flow surface, when an abrupt change in 
flow direction is indicated by the structural characteristics values of the latest 
grouped cell, etc. 

Termination Identification 

The points at which the geologic strata terminate or become relatively 
close to each other may also be identified. As discussed above, identifying 
termination points is an important activity associated with the location and 
evaluation of stratigraphic traps. If the dip and azimuth of a particular 
geologic stratum is traced and we detect where this stratum terminates 
against another stratum, an onlap, downlap, toplap, or erosional truncation is 
likely and these features are associated with potential stratigraphic 
hydrocarbon traps. 

If two flow surfaces intersect, this may correspond to a geologic strata 
termination interface. However, due to interference in the seismic signal and 
due to the smoothing in the dip/azimuth estimate, flow surfaces tend to 
converge rather than intersect at the termination interfaces. One way to 
alleviate this is to assign a width or window to each flow surface. If a flow 
surface intersects an area covered by the window around another flow 
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surface, the intersection is interpreted as a termination. Terminations may be 
detected practically by producing a binary volume with the flow surfaces while 
the flow surfaces are generated. If a flow surface approaches a cell marked 
as visited by an earlier flow surface, this cell is marked as a geologic stratum 
termination cell. 

It is also noted that terminations have proven to be good indicators of 
other geologic features and seismic signal characteristics than terminations, 
e.g. faults. The identified termination cells may therefore be used for 
purposes other than solely for stratigraphic interpretation. 

Stratigraphic Fades Identification 

The structural characteristics values may also be used in the 
identification of the stratigraphic fades, i.e. the internal structure or texture of 
the strata in the subsurface area. Several alternative processes may be used 
to compute and represent this texture information. Some will yield an 
attribute value for a confined sub-volume, while others will yield one attribute 
value for each seismic sample. The dip/azimuth estimation technique and the 
flow surface descriptions discussed above provide powerful primitives for 
generating fades texture attributes. 

Attributes can be created, for instance, that highlight chaotic regions of 
the seismic signal. One such method uses the smoothed magnitude of the 
gradient in the dip and azimuth estimates, normalized with respect to the 
smoothed magnitude of the seismic. If the dip/azimuth estimate is reliable, 
this measure will yield a high value. Chaotic regions will typically not yield 
reliable estimates. An alternative is the reliability measure of the dip and 
azimuth estimate discussed above. Similar attributes can be calculated by 
measuring the regularity of the flow lines. For each sample position, two flow 
lines can be generated starting from a distance above and below this sample 
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position. The attribute value at this sample position is then set to be the 
standard deviation of the change in distance between the flow lines. In 
regions with chaotic signals, the distance between the flow lines will change a 
great deal, and the standard deviation will be high. 

Attributes can also be developed that highlight other types of 
statigraphic fades patterns, such as divergent, straight, and wavy parallel 
patterns. One such type of attribute can be calculated by initiating flow lines 
of a certain length above and below each sample of the seismic data. Then 
the distance between the flow lines' respective left and right end points, d, and 
d r can be calculated. Attribute values can then be calculated as \d, -rf r |. 
Divergent signals will lead to flow lines with consistent differences in the 
differences in the distance between the flow line end points, whereas parallel 
signal configurations will lead to fairly constant distances. A similar method 
computes the attribute value as the mean value of the flow lines' separation 
gradient normalized by the standard deviation of the separation gradient. 
Divergent flow lines will have a somewhat constant gradient with nonzero 
mean and small standard deviation. 

It is easy to see how these techniques can yield a good separation of 
the parallel from the divergent signal characteristics. Used in combination 
with each other or with the chaos attributes its is then easy to discriminate 
between divergent, parallel and chaotic stratigraphic fades. 

Represent/Display Results 

The cells grouped in the Group Cells step 20 collectively provide a 
granular representation of the identified geologic feature. While this type of 
representation will be sufficient for many purposes, the identified subsurface 
geologic feature may show undesired characteristics, such as holes or 
inaccurately grouped cells, and several different subsurface features may be 
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connected into a single composite feature in the data volume. This latter 
problem is particularly likely if the feature being identified is a fault and the 
seismic data being interpreted has two or more faults that cross each other. 

These issues may be addressed in the Represent/Display Results step 
22. In this step, the grouped cells may be clustered or subdivided and the 
grouped cells or these clustered portions of the grouped cells may be 
represented as vectorized surfaces in three dimensions, such as a plane: 

an, + bn 2 + cn 3 = d, 

where n 1t n 2i and n 3 represent the three axis coordinates, or a higher order 
polynomial surface, such as 

an 2 , + bn, + cn 2 + dn 3 3 + en 3 = f. 

If the surface shape is too complex for such a model, each surface 
may be subdivided into multiple connected patches of polynomial models. 
The parameters of this surface representation may be estimated using least 
squared error approximation or other parameter estimation techniques known 
to those skilled in the art. 

If less than all of the grouped cells are to be represented, the cells 
must be suitably clustered, i.e. which samples correspond to which surfaces 
must be determined. In some cases, the clustering and vectorizing processes 
may be performed as a sequential process. In other cases, better results 
may be obtained using an iterative method where the results of the 
vectorizing process are used to refine the results of the clustering process. 

Several alternative methods for clustering the grouped cells may be 
used in the inventive method. In one method, the grouped cells are 
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subdivided and connected components are extracted from these subdivided 
units. Connected components are sets of cells that are connected within a 
geometrical neighborhood. A parametric surface representation of these 
connected components is made and the degree of fit of this parametric 
surface representation with respect to all of the samples in the larger volume 
is measured. Next, the cells with best degree of fit are used to recalculate the 
surface parameters, then new degrees of fit are computed, etc. The last part 
of this process is iteratively repeated until a sufficient degree of fit is obtained. 

Other clustering schemes may also be employed including letting an 
interpreter select (portions) of the connected components to start from; 
estimating parameters, and iteratively fitting computations and reestimations 
as above. Detailed information on general clustering techniques may be 
found in Pattern Recognition: Statistical, Structural and Neural Approaches by 
Robert Schalkoff, John Wiley, 1992, incorporated herein by reference. 

The grouped, clustered, and/or vectorized cells may then be suitably 
displayed as shown in Figures 4 through 9. Figures 4 and 5 display the types 
of subsurface fault features that may be identified by the present method. In 
Figure 4, section 50 is a sectional view through a seismic data cube with fault 
features 52 showing the types of faults that may be identified by the present 
method. Figure 5 shows a perspective view of the entire seismic data cube 
60 as well as perspective view of faults 62 identified by the present method. 
In the processing used to create Figures 4 and 5, the parameters used were 
<7, =3.0, o 2 =3.0, <r 3 =1.0, <x 4 =6.0 and a 5 =5.0. 

In this example, the cells were grouped into blocks of size 40x40x40 
when producing the vectorized fault representations shown in Figures 4 and 5 
and the parametric surface representation used was: 

p t n 3 \ + p 7 n \ + /? 3 /7, + p A n 7 + p % n y i + p t n 7 i + /? ? w 3 = 1 , 



PCT/IB99/01040 

34 

where the p/s are the surface parameters. These parameters were estimated 
by least squares approximation. 

Figures 6 and 7 show the types of seismic reflector features that may 
be identified by the present method. In Figure 6 f section 70 is a sectional 
view through a seismic data cube (a collection of seismic data obtained from 
a particular subsurface area). First seismic reflector feature 72 and second 
seismic reflector feature 74 represent a pair of seismic reflectors of the type 
that may be identified by the present method. Derrick 75 shows a 
hypothetical well location and first wellbore location 76 and second wellbore 
location 78 represent a pair of locations within the well that could have acted 
as seed points for the location of first seismic reflector feature 72 and second 
seismic reflector feature 74. Figure 7 shows the types of three-dimensional 
seismic reflector features that may be identified by the present method. The 
first reflector surface 80 and second reflector surface 82 shown in Figure 7 
may represent a more complete 3D view of the seismic reflectors shown in 
cross section in Figure 6. 

The type of seismic reflector termination information that may be 
determined using the inventive method is shown in the 2D cross section in 
Figure 8. In the processing used to create Figure 8, flow lines were initiated 
at every row in every 10th column, using flow lines having a width of three cell 
or sample positions. If a flow line intersected the window around another flow 
line which was initiated at least five rows away on the same column, a 
termination density map was updated by adding one to that cell. The 
momentum concept was also used, with the last 31 angular estimates along 
the flow line being averaged. 

Figure 9 is a cross-sectional view through a fades texture attribute 
cube prepared in accordance with the inventive method that highlights chaotic 
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signal patterns as discussed above. Section 90 identifies a particular 
subsurface area 92 where the seismic data signal is highly chaotic. 

The present method is particularly suited for autonomous or semi- 
autonomous seismic data interpretation, where steps 18, 20, 22, and 24 are 
performed with little or no manual operator input. 

It will be readily understood that the steps and processes associated 
with the disclosed embodiment of the present method are capable of a wide 
variety of alternative implementation methods. The entire seismic data cube 
may be examined or only a particular area within the cube. Many of the steps 
may be performed iteratively. 2D methods may be used to identify 2D 
portions of the subsurface geologic feature and these 2D portions may then 
be aggregated and interpolated to produce a 3D representation of the 
subsurface geologic feature. The present method is in no way limited or 
restricted to the particular order of steps described in the preferred 
embodiment above. 
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CLAIMS; 

1 . A method of identifying a subsurface geologic feature using seismic 
data obtained from a subsurface area, said method comprising the steps of: 

subdividing said subsurface area into a plurality of discrete cells; 
each cell representing a different region within said subsurface area; 

associating said cells with portions of said seismic data that 
image subsurface regions represented by said cells; 

determining structural characteristics for a plurality of said cells 
using said seismic data; and 

grouping together nearby said cells based on said structural 
characteristics to identify said subsurface geologic feature. 

2. A method according to Claim 1 , wherein said subsurface geologic 
feature comprises a fault. 

3. A method according to Claim 1, wherein said subsurface geologic 
feature comprises a seismic reflector. 

4. A method according to Claim 1 , wherein said subsurface geologic 
feature comprises a hydrocarbon reservoir boundary. 

5. A method according to Claim 4, wherein said hydrocarbon reservoir 
boundary comprises a hydrocarbon/water interface. 
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6. A method according to Claim 1, wherein said grouping step comprises 
selecting nearby cells having essentially horizontal local seismic reflector 
inclinations. 

7. A method according to Claim 1 , wherein said determined structural 
characteristics comprise an estimate reliability measure. 

8. A method according to Claim 7, wherein said grouping step comprises 
grouping together nearby cells having relatively high estimate reliability 
measures. 

9. A method according to Claim 7, wherein said grouping step comprises 
grouping together nearby cells having relatively low estimate reliability 
measures. 

10. A method according to Claim 1 , wherein said grouping step comprises 
the step of selecting a seed cell from said plurality of cells, said seed cell 
representing a point on a first seismic reflector. 

11. A method according to Claim 10, wherein said grouping step further 
comprises the step of selecting additional cells at increasing distances from 
said seed cell using said local seismic reflector inclinations, said additional 
cells representing additional points on said first seismic reflector. 

12. A method according to Claim 1 1 , wherein said selected cells when 
aggregated with neighboring selected cells produce localized seismic reflector 
areas generally in alignment with said local seismic reflector inclination 
estimates of said cells through which said seismic reflector passes. 

13. A method according to Claim 10, further including the step of extracting 
a flow surface using said seed cell and said structural characteristics. 
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14. A method according to Claim 13, wherein said flow surface is extracted 
using momentum. 

15. A method according to Claim 13, wherein a plurality of said flow 
surfaces are extracted. 

16. A method according to Claim 15, wherein cells where said extracted 
flow surfaces become relatively close together are identified as seismic 
reflector terminations. 

17. A method according to Claim 1, wherein said seismic data has been 
time domain to depth domain converted prior to said determining step. 

18. A method according to Claim 1, wherein said seismic data has been 
decomposed and reconstructed prior to said determining step. 

19. A method according to Claim 1, wherein said cells are thinned in said 
grouping step. 

20. A method according to Claim 1, wherein said grouped together cells 
are subdivided using a clustering process. 

21 . A method according to Claim 1 , wherein a plurality of said grouped 
cells are represented as a vectorized surface. 

22^ A method according to Claim 1, further comprising the step of 
displaying said grouped cells. 

23. A method according to Claim 1, wherein said seismic data associated 
with said cells are associated with a single subsurface imaging point. 
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24. A method according to Claim 1, wherein said determining step 
comprises orientation selective filtering said seismic data. 

25. A method according to Claim 24, wherein said orientation selective 
filtering comprises estimating gradients of local seismic reflectors. 

26. A method according to Claim 25, wherein said gradient estimation is 
performed by filtering said seismic data with a derivative of Gaussian filter. 

27. A method according to Claim 26, wherein said gradient estimates are 
smoothed. 

28. A method according to Claim 27, wherein said smoothing is performed 
by filtering said gradient estimates with a Gaussian low-pass filter. 

29. A method according to Claim 1 , wherein said structural characteristics 
comprise a fault characteristic that quantifies a likelihood that a fault passes 
through said cell. 

30. A method according to Claim 29, wherein said grouping step 
comprises grouping together nearby cells having relatively large fault 
characteristic values. 



31 . A method according to Claim 1 , wherein said cells are three- 
dimensional volume elements. 

32. A method according to Claim 31 , wherein said structural characteristics 
comprise dip and azimuth values. 
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33. A method according to Claim 32, wherein said dip and azimuth values 
are determined by principal component analysis of smoothed gradient 
estimates. 



34. A method according to Claim 1, wherein said cells are two-dimensional 
area elements. 

35. A method according to Claim 34, wherein said structural characteristics 
comprise plunge angles. 

36. A method according to Claim 1, wherein said structural characteristics 
comprise stratigraphic facies attributes. 

37. A general purpose computer programmed to implement the method of 
any one of Claims 1 to 36. 

38. A computer program on a carrier containing instructions for causing a 
general purpose computer to implement the method of any one of claims 1 to 
36. 
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